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In this Brief Report, we present an algorithm for calculating the elastic Lattice Greens Function 
of a regular lattice, in which defects are created by removing lattice points. The method is 
computationally efficient, since the required matrix operations are on matrices that scale with 
the size of the defect subspace, and not with the size of the full lattice. This method allows the 
treatment of force fields with multi-atom interactions. 



Lattice Green's function methods provide a very com- 
putationally efficient way of handling the long-ranged 
stress fields around defects in materials. The purpose 
of this Report is to generalize pthe lattice Green's func- 
tion method of Thomson et alH to problems in which a 
set of lattice points has been completely removed from 
the problem. Such problems may be the modeling of 
the elastic fields in crystalline systems in which cavities 
are created. An example of such systajBS that have been 
modeled in this way are blunt cracksoS. This opens for 
the application of Green's function techniques to a np 
class of problems: The previously published methoc 
has been used to study lattice defects whexe^fonnections 
between lattice points have been alterectl'B~Q (for exam- 
ple leading broken to introduce cracks). The present 
extension of the method permits the study of problems 
where lattice points (atoms) have been removed from the 
system. The methodology was used in the calculations of 
Refs. H and ^, but the derivation has not been presented 
in the literature yet. 

In the lattice Green's function method, one begins with 
a force constant matrix defined by 
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where E is the total energy, is the displacement in 
the i direction of atom a and summation over repeated 
indices is implied, is thus the negative of the force 
in the i direction on atom a when atom b is displaced a 
unit distance in the direction j: 
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For a simple pair potential with radial forces between 



atoms, ^ is given by 
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where /(r) is the force between the two atoms, r is the 
vector between the equilibrium positions of atoms a and 
b and e.; is a unit vector in the i direction. The rf''/ \r°-''\ 
are a projection operators. 

The Greens function gives the formal solution of equa- 
tion (0) for u: 



G = * 



(4) 
(5) 



It is inconvenient to label the degrees of freedom by two 
indices (the atom and the component of the displace- 
ment). Therefore, we introduce a single labeling of all 
degrees of freedom, leading naturally to a matrix nota- 
tion where a component of G, G^, is the response by 
the degree of freedom i to a force acting on the degree of 
freedom j. _ 

Thomson et al. demonstratetJ how the Greens function 
of a regular lattice can be calculated in Fourier space in 
a time that scales linearly with the number of lattice 
points {Niat)- The Fourier-space approach breaks down 
when defects are present, breaking translational symme- 
try. Thomson et al. showed that defects can be treated 
without a complete "brute- force" inversion of equation 
(^ , by the use of multiple-scattering theory. If the val- 
ues of only a small number of "bonds" (i.e. elements in 
the $ matrix) are changed, the new Greens function G 
can be calculated from the perfect-lattice Greens function 
G*^ using a Dyson equation: 



G = G"(l - 5*G" 



(6) 



where 5^ is the change in the force constant matrix. This 
calculation scales as N^ef^ where Ndef <C Niatuce is the 
number of rows or columns in (5$ containing non-zero 
elements, i.e. the number of degrees of freedom of atoms 
where one or more bonds are modified, and Niatuce is the 
total number of degrees of freedom in the lattice. 

When atoms are removed from the system (see fig- 
ure |l|) , it is important to be able to calculate the Green's 
function as efficiently as when only bonds are changed, 
i.e. the Greens function should be calculated starting 



1 




FIG. 1. A hexagonal lattice, where 12 atoms have been 
removed. The removed atoms and their bonds are shown 
with dashed lines. Border atoms that have lost one or more 
neighbors are shaded. The corresponding diagonal elements of 
5^ are non-zero, since the on-site force constant has changed: 
when such an atom is displaced it "feels" a different force since 
there is no force from the missing neighbors. In the text the 
following classification scheme is used: The removed atoms 
are in class B, the remaining atoms are in class A. Class A 
is further subdivided into two classes: the grey atoms, having 
lost one or more neighbors, are in class Ai, and the white 
atoms are in class A2. 

from the perfect-lattice Greens function G°, and not di- 
rectly from equation Let us divide the degrees of 
freedom in two classes: A are the degrees of freedom 
of atoms that are "kept" , and B are those that are re- 
moved from the system. M^a is then the sub-matrix 
of M obtained by keeping only columns and rows corre- 
sponding to degrees of freedom in class A. In a similar 
ways the other sub-matrices M-bb, Mas, and M.ba are 
defined. Further, let G* be the Green's function of the 
new system. G* cannot be calculated by breaking all 
bonds to the atoms to be removed, and then applying 
the Dyson equation (|^), since the resulting Greens func- 
tion is singularll3. 

The force-constant matrix entering G* is perturbed 
relative to the full force-constant matrix in two ways. 
First, it has smaller dimension in the sense that it only 
contains atoms in the ^-region. Second, the on-site force 
constants (corresponding to the force on an atom result- 
ing from its own motion) of the A-region border atoms 
near the cavity are changed, since these are determined 
by the neighboring atoms (terms in the sum in equation 
(pq) are missing). For example, if an atom has no neigh- 
bors, then it has a vanishing on-site force constant. Then 
G* is given by 

G* = ($AA + S^aaV' , (7) 

where 6^aa corresponds to the change in the on-site 
force constants of the border atoms. For clarity, we de- 
fine the auxiliary Green's functions G'^ = €>^^ for the 
perfect crystal, and G' = ($aa)^^ for the crystal with 
atoms removed but no changes in the force constants 
of the remaining atoms. (Note that G' — ($aa)^^ 



(*-^)aa = Gaa)- The calculat ion of G* then goes in 
two steps. We first calculate G', as an intermediate step: 

G' = (*aa)-^=(((GO)-^)^J ' (8) 

— ^AA ^ ^abK^bb) ^ba 

(The last equality is proven in the appendix). We then 
calculate G* with a Dyson equation: 

G*(*aa + '^*aa) = 1 (9) 

which gives 

G* = G' (1 - G'^$aa) ~^ (10) 

Equations ^ and ([lo| ) are the central results of this 
Report. We note that they can be used to describe sys- 
tems with many-atom interactions such as the Eptbed- 
ded Atom MfithodB, the Effective Medium TheoryEl, the 
glue modeltj etc. In these methods, the force constants 
between different atoms in the A region also change be- 
cause the effective interatomic potentials are environmen- 
tally dependent. This effect can be included straightfor- 
wardly by modification of the 6^aa term. 

Finally, let us show that this method is computation- 
ally feasible. Let the degrees of freedom that are kept be 
separated into two classes, Ai and A2, where Ai is the 
set of degrees of freedom of atoms where the correspond- 
ing elements of 6^ is non-zero, or where we are going to 
need the Greens function; the latter would include atoms 
on which nonlinear forces or loading forces will eventu- 
ally be placed. Let n{S) be the number of degrees of 
freedom of the atoms in set S. Typically, n{Ai) will be 
a few hundreds or Jess, whereas A2 contains millions of 
degrees of freedomBB. n{B) is also typically a few hun- 
dreds. Equation (||) then consists of the inversion of a 
n{B) X n{B) matrix, and two n{Ai) x n(B) matrix mul- 
tiplication. Since S^A2A2 — equation (|lO| ) consists of 
a n{Ai) X n(Ai) matrix inversion and a similar multipli- 
cation. With the relatively small sizes of Ai and B, this 
is clearly computationally practical, whereas the direct 
matrix inversion in equation (^ is not. See also the dis- 
cussion of the computational burden of equation ([l^) in 
Thomson et alE if there are many atoms in Ai where 
is zero. 

In conclusion, we have shown how to generate Lattice 
Green's functions for lattices where the defects are not 
limited to perturbations of a small number of bonds, but 
where lattice points have been removed. The Green's 
function is generated in a computationally efficient way, 
and has already beea used to simulate the elastic fields 
near a blunted crackEl. The method is not limited to elas- 
tic fields in crystal lattices, but to all problems where a 
lattice Green's function may be useful, and where a small 
number of lattice points are removed from the problem. 
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APPENDIX 



We here prove that 
1 
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= Maa-Mab{Mbb) Mba (11) 



Proof: Let N = M"^ We then have MN = 1, which 
can be split into four parts: 



Maa^aa + Mab^ba = 1, 
Maa^ab + Mab'Nbb - 0, 
Mba'!^aa + MbbNba = 0, 
MbbNbb + MssNsij = 1. 



From (12c) we get 

Nba = -(Mbb)"^MbaNaa. 



(12a) 
(12b) 
(12c) 
(12d) 



(13) 



When inserting in equation (12a) and multiplying from 
the right with (N^a), we get 



Maa - Mab{Mbb) ^Mba = (Naa)" 



(14) 



and since (N^a) ^ ((M ^)aa) ^, equation ( p]) has 
been proven. 
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